Download US Census data on natural resource dependent jobs (total # and median income) within counties and summarize across NEP boundaries:
nep_sf <- st_read(dsn = "./data-raw", layer = "NEP_Boundaries10162018")
## Reading layer `NEP_Boundaries10162018' from data source `D:\rprojects\anep_bea_stats\data-raw' using driver `ESRI Shapefile'
## Simple feature collection with 28 features and 8 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -124.7411 ymin: 18.31383 xmax: -65.87813 ymax: 48.99944
## epsg (SRID): 4269
## proj4string: +proj=longlat +datum=NAD83 +no_defs
states <- c("AL", "AZ", "AR", "CA", "CO", "CT", "DE", "FL", "GA",
"ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD",
"MA", "MI", "MN", "MS", "MO", "MT", "NE", "NV", "NH", "NJ",
"NM", "NY", "NC", "ND", "OH", "OK", "OR", "PA", "RI", "SC",
"SD", "TN", "TX", "UT", "VT", "VA", "WA", "WV", "WI", "WY",
"PR")
metrics <- load_variables(2017, "acs5", cache = TRUE) %>%
filter(grepl("Fishing|fishing", label))
totaljobs_sf <- get_acs(geography = "county", variables = c(Ag_For_Fish_Hunt_Min = "C24050_002"),
state = states, geometry = TRUE)
nr_medinc_sf <- get_acs(geography = "county", variables = c(Med_Income = "B24031_003"),
state = states, geometry = TRUE)
ustotaljobs_sf <- st_transform(totaljobs_sf, st_crs(nep_sf))
usnrmedinc_sf <- st_transform(nr_medinc_sf, st_crs(nep_sf))
nep_jobs_intersects <- st_intersects(nep_sf, ustotaljobs_sf)
nep_medinc_intersects <- st_intersects(nep_sf, nr_medinc_sf)
nep_sel_sf <- ustotaljobs_sf[unlist(nep_jobs_intersects),]
nep_sel2_sf <- nr_medinc_sf[unlist(nep_medinc_intersects),]
nep_jobs <- st_join(nep_sf, nep_sel_sf, join = st_intersects) %>%
sf::st_buffer(dist = 0)
nep_nrmedinc <- st_join(nep_sf, nep_sel2_sf, join = st_intersects) %>%
sf::st_buffer(dist = 0)
nep_jobs_sum <- nep_jobs %>%
select(NEP_NAME, estimate) %>%
group_by(NEP_NAME) %>%
summarise(jobs = sum(estimate))
nep_medinc_sum <- nep_nrmedinc %>%
select(NEP_NAME, estimate) %>%
group_by(NEP_NAME) %>%
summarise(medinc = median(estimate))
Plot US Census natural resource dependent jobs data within the NEP boundaries:
pal <- colorNumeric(palette = "viridis", domain = nep_jobs_sum$jobs, n = 10)
nep_jobs_sum %>%
st_transform(crs = "+init=epsg:4326") %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addPolygons(popup = ~ paste(NEP_NAME,"</br>","Jobs = ",jobs),
stroke = FALSE,
smoothFactor = 0,
fillOpacity = 0.7,
color = ~ pal(jobs)) %>%
addLegend("bottomright",
pal = pal,
values = ~ jobs,
title = "Natural Resource Dependent Jobs</br>(Ag.,Forest.,Fishing,Hunting,Mining)",
opacity = 1)
Plot US Census natural resource dependent jobs, median income data within the NEP boundaries:
pal2 <- colorNumeric(palette = "viridis", domain = nep_medinc_sum$medinc, n = 10)
nep_medinc_sum %>%
st_transform(crs = "+init=epsg:4326") %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addPolygons(popup = ~ paste(NEP_NAME,"</br>","Median Income = ",medinc),
stroke = FALSE,
smoothFactor = 0,
fillOpacity = 0.7,
color = ~ pal2(medinc)) %>%
addLegend("bottomright",
pal = pal2,
values = ~ medinc,
title = "Natural Resource Dependent Jobs</br>(Median Annual Income)",
opacity = 1)
Import US BEA data on natural resource dependent jobs (total county income) and summarize across NEP boundaries:
nr_income <- read.csv("./data-raw/CAINC5N__ALL_AREAS_2001_2017.csv", header = TRUE) %>%
filter(LineCode==100,str_detect(GeoFIPS, "....0")==FALSE) %>%
mutate(GEOID = as.character(GeoFIPS),
YR2017 = as.character(X2017)) %>%
select(GeoName, GEOID, YR2017)
nr_income$GEOID <- gsub(" ", "", nr_income$GEOID, fixed = TRUE)
nr_income$Y2017 <- as.numeric(nr_income$YR2017)
nr_inc_sf <- left_join(nep_sel_sf, nr_income, by = c('GEOID' = 'GEOID'))
nep_totinc <- st_join(nep_sf, nr_inc_sf, join = st_intersects) %>%
sf::st_buffer(dist = 0)
nep_inc_sum <- nep_totinc %>%
select(NEP_NAME, Y2017) %>%
group_by(NEP_NAME) %>%
summarise(totinc = sum(Y2017, na.rm = TRUE))
nep_inc_sum$totinc <- ifelse(nep_inc_sum$totinc==0,NA,nep_inc_sum$totinc*1000)
Plot US BEA natural resource dependent jobs total income estimates within the NEP boundaries (where available):
pal3 <- colorNumeric(palette = "viridis", domain = nep_inc_sum$totinc, n = 10)
nep_inc_sum %>%
st_transform(crs = "+init=epsg:4326") %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addPolygons(popup = ~ paste(NEP_NAME,"</br>","Total Income (2017) = ",totinc),
stroke = FALSE,
smoothFactor = 0,
fillOpacity = 0.7,
color = ~ pal3(totinc)) %>%
addLegend("bottomright",
pal = pal3,
values = ~ totinc,
title = "Natural Resource Dependent Jobs</br>(Total Annual Income, if reported)",
opacity = 1)